function [LE,MU]=Dynkw(B,nf,ny,lpd,bcrit,track,schur,nb)

% This program computes the LE,MU matrices that are used
% to solve the "fundamendal dynamic equation" of the model,
% which has been isolated in the program redkw.m
%
% This equation takes the form
%
%    Ed(t+1)|I(t) = W d(t) + psid(F) Ex(t)|I(t)
%
% where d(t) is a vector of dynamic variables in the reduced
% system.  The program redkw.m found "nf" elements of the y(t)
% vector that are excluded from this dynamic system.  Thus, there 
% are nd=ny=nf elements of d(t).  
%
% To solve this system, we must find a "diagonalization" of it, 
% as in Blanchard and Kahn, that allows for the decoupling of the
% dynamic variables that are unstable from those that are stable.
% There are two methods allowed: an eigenvector method (with or 
% without the balancing of the eigenvalues) and a schur decomposition
% method.  The latter is numerically most stable and accurate, but
% may be slower.
if nargin<7;
 schur=0;
end

% dynkw.m
% solution for the LE,MU matrices to be used in dynamic solution
% programs (relies on function lus.m)

v=version;

if (eval(v(1))>=5)
   np=length(lpd); 
else
   np=max(size(lpd));
end

if (nargin<8)
   nb=0;
   if (track==1)&(schur==0)
      disp(' ')
      disp('dynkw.m: default is balanced eigenvalues')
      disp('set arg 8 to nb=1 for nonbalanced eigenvalues')
   end
end

% Selection of the W matrix as the last nd=ny-nf rows and columns of B

nd=ny-nf;   %nd= total number of state variables (predetermined and nonpredetermined)

W=B(nf+1:ny,nf+1:ny);

if (schur~=1)
   
% The subroutine eigl.m computes the left eigenvectors of W, using
% a balance (nb=0) or no balance (nb=1) option as indicated by the 
% researcher in the call to dynkw.m


[LE,MU]=eigl(W,nb);



else

% The schurl.m program yields a decomposition, Q'*W*Q=T
%     Where T is LOWER triangular  (Note lower and not upper)
%     and the diagonal entries of T are the ordered
%     eigenvalues of  W with T[i,i] >= T[j,j], for i < j.


[Q,T]=schurl(W);
LE=Q';
MU=T;

end

% In either case, we now have matrices that isolate the unstable
% eigenvalues first and then isolate the stable eigenvalues.
% Determine the number of eigenvalues greater than critical level

nus=sum(abs(bcrit*diag(MU))>1);

% Select the portion of the LE matrix that isolates these eigenvalues.

if (nus>0)
 LE=LE(1:nus,:);         
 MU=MU(1:nus,1:nus);
else
 LE=[];
 MU=[];
end

% Tests of the decomposition

% Test 1: are there the same number of unstable roots
%         as there are nonpredetermined variables in 
%         the reduced system?

if (track==1)
   disp(' ')
   disp(['dynkw.m: critical value for instability is ', num2str(1/bcrit)])
   disp(['which corresponds to bcrit = ',num2str(bcrit)])
   disp(' ')
end

if (nus>nd-np)
  disp(' ')
  disp(['Using a critical value of bcrit =',num2str(bcrit)])
  disp('There are too many unstable roots, i.e., roots such that')
  disp('bcrit*mu >=1.  Hence, there is no stable solution to the model')
  disp('This can be due to an inappropriate choice of bcrit, an intrinsically')
  disp('unsolvable model, or an error in the specified system')
  disp('Use debugkw.m to look into system errors, which could include')
  disp('the wrong number of predetermined variables')
  disp('  ')
  disp(['Dimension of reduced system            ', num2str(np)]) 
  disp(['Number of nonpredetermined variables:  ', num2str(nd-np)])
  disp(['Number of unstable roots               ', num2str(nus)])
  disp('  ')
  disp('All results that may be obtained after this point are meaningless')
  pause
end

if (nus<nd-np)
  disp(' ')
  disp(['Using a critical value of bcrit =', num2str(bcrit)])
  disp('There are too few unstable roots, i.e., roots such that')
  disp('bcrit*mu >=1. Equivalently, there are too many stable roots')
  disp('Hence, there are multiple stable RE solutions to the model')
  disp('This can be due to an inappropriate choice of bcrit, a model that')
  disp('cannot be uniquely solved, or an error in the specified system')
  disp('Use debugkw.m to look into system errors, which could include')
  disp('the wrong number of predetermined variables')
  disp('  ')
  disp(['Dimension of reduced system            ', num2str(np)]) 
  disp(['Number of nonpredetermined variables:  ', num2str(nd-np)])
  disp(['Number of unstable roots               ', num2str(nus)])
  disp('  ')
  disp('All results that may be obtained after this point are meaningless')
  pause
end

if (nus>0) % If there are unstable roots
   
% Test 2: is it possible to associate the unstable canonical
%         variables and the nonpredetermined variables?

% Note: here you'd want the number of nonpredetermined state variables to be equal to the number of unstable roots




Llam=LE(:,1:nus);
rLlam=rank(Llam);

if (nus>rLlam)
 disp(' ')
 disp('dynkw.m: rank condition on Llam is violated')
 disp('number of unstable roots')
 disp(nus)
 disp('rank of Llam')
 disp(rLlam)
 disp('Llam')
 disp(Llam)
 disp('model cannot be solved: program should be terminated')
 disp('run debugkw.m to check system') 
 pause
end

end